Numerical solution for circular tunnel excavated in strain-softening rock masses considering damaged zone

Despite the extensive investigation on the stress and displacement distributions in tunnels, few have considered the influences of the damaged zone around a tunnel on the strength and stiffness parameters of the surrounding rock, including the gradual variation in the damaged factor D and dimensionless damaged radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho^{{\text{d}}}$$\end{document}ρd, under the effect of excavation disturbance. In this paper, a numerical solution is presented for the stresses and displacement of a circular tunnel excavated in a Hoek–Brown rock mass considering the progressive destruction of the damaged zone. First, the results obtained using the proposed algorithm are compared with those obtained using other numerical solutions, demonstrating a high degree of accuracy through some examples. Second, the influences of the damaged factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D$$\end{document}D and dimensionless damaged radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho^{d}$$\end{document}ρd on the stresses, radial displacement, plastic radii, and ground response curve are investigated. The results show that, as the damaged factor D increases, the radial displacement and plastic radius increase, whereas the tangential stress decreases. Both the plastic radius and displacement decrease with decreasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho^{{\text{d}}}$$\end{document}ρd. This shows that the damaged factor D has a significant effect on tunnel convergence.


Problem definition and methodology
shows a schematic of a circular tunnel of radius b excavated in homogeneous and isotropic rock masses under the effect of a far-field hydrostatic stress σ 0 . In Fig. 1, R p and R d represent the plastic and damage radii; σ θ and σ r represent the tangential and radial stresses around the tunnel, respectively. In such cases, with a decrease in the internal support pressure p i , an elastic deformation of the surrounding rock can occur. When the support pressure p i is less than the critical pressure p ic , a plastic zone may be formed around the tunnel opening 22 . With the formation of the plastic zone, an excavation damage zone with a radius R d will gradually form inside the plastic zone 25,30 , and the damaged factor D gradually increases from the depth of the damage zone, reaching maximum at the critical point in the residual zone. The list of symbols used are shown in the following Table 1. www.nature.com/scientificreports/ Yield criterion and strain-softening behavior of materials. The strength-weakening behavior can be modeled according to the theory of plastic mechanics, which can help derive the elastoplastic deformation process of the surrounding rock 31 . Based on this theory, both the failure criterion F and the plastic potential function G depend on the stress state and the softening parameter η p of the rock mass. Therefore, it is important to select appropriate softening parameters and yield criteria 32 . In this study, the plastic shear strain η p was used as the deviator plastic strain parameter 18 : where ε p θ and ε p r are the tangential and radial strains representing the major and minor principal plastic strains ε p 1 and ε p 3 , respectively. Therefore, Eq. (1) can be rewritten as follows: The yield criterion can be expressed as follows: The nonlinear H-B instability criterion is similar to that in the case of rock plastic deformation. Therefore, H(σ r , η p ) can be expressed as 30 : where σ c represents the uniaxial compressive strength of the rock mass; m , s , and α represent the strength parameters of the H-B rock mass, which can be obtained using the geological strength index (GSI) 33 . Therefore, the peak and residual values of the strength parameters can be expressed as 26 : www.nature.com/scientificreports/ where the subscripts "p" and "r" denote the peak and residual values of the surrounding rock parameters m, s, α and GSI, respectively. In addition, m i and D represent the material constant of the rock under intact conditions and the degree of rock mass damage, respectively. The residual GSI r of the rock mass can be calculated using the following formula 34 .
In addition to estimating the GSI of the undamaged and damaged rock masses, estimating the deformation modulus of the rock mass is an important part of the deformation calculation. Therefore, the displacement analysis of the tunnel also requires estimating the deformation modulus of the excavated rock mass. Many scholars have shown that the deformations modulus is not a constant 35 . Based on the database of rock mass deformation modulus measurements, the following expression has been proposed for estimating the rock mass modulus 4 : where E 0 is the intact rock deformation modulus.
Plastic potential function of the material. Selecting an appropriate plastic potential function has an important influence on the calculation of the plastic strain process. In geo-mechanics, the Mohr-Coulomb type of plastic potential function is widely applied and implemented, which can be written as 25 : where k η p is the coefficient of dilation and is computed using Eq. (14).
where ϕ is the expansion angle of the material 17 . In this study, the increment relationship between the radial and tangential plastic strains can be obtained according to the non-associated flow rule 21 .
Evolution of material parameters in different zones. Based on field experience, a new elastoplastic damage piecewise curve strain softening behavior model is established for undamaged and damaged rock masses, as shown in Fig. 2. The stress and strain behavior follows the piecewise linear softening behavior in the plastic undamaged and damaged zones, which are represented by solid red and blue lines in Fig. 2, respectively.
As mentioned above, the strength parameters of the material vary depending on the softening functions employed for the different zones. The evolution of the strength parameter of the rock mass can be described using the piecewise functions of the deviatoric plastic strain η p in the plastic undamaged and damaged zones, and calculated using Eqs. (16) and (17), respectively.
where ω u (η) represents one of the rock mass parameters,m(η) , s(η) , α(η) , GSI(η) and ϕ(η) in the undamaged zone and ω d (η d ) represents one of the rock mass parameters, m d (η) , s d (η) , α d (η) , GSI d (η) , D(η) and ϕ d (η) in the damaged zone 25 . In the above formula, D p = 0 represents the boundary at the beginning of the damaged zone, and D r represents the maximum damage degree of the damage zone. η * is the critical plastic shear strain representing the starting point of the residual behavior and needs to be determined experimentally 32 . The subscripts "p" and "r" indicate the peak and residual values 21 , respectively, and the superscripts "u" and "d" represent the undamaged and damaged zones, respectively.

Basic numerical formulations for strain-softening model
Preliminaries. In this section, the complex stress-strain relationship is derived based on the newly defined medium. Evidently, the radial stress σ r = σ R at r = R p (outer boundary of the plastic zone) in Fig. 1. Assuming that the plastic zone is composed of n concentric annuli, as shown in Fig. 3, the ith annulus is surrounded by two circles of normalized radii At the elastoplastic boundary, where i = 1, according to the theory of elasticity, the stress and strain can be obtained as follows 14 : where σ and ε represent the stress and strain, respectively, and the subscripts "r" and " θ " represent the radial and circumferential directions, respectively. E (1) is the Young's modulus, and v is the Poisson's ratio of the intact rock mass at the elastoplastic boundary 3 .
Based on previous research 21,32,36,37 , the radial stress increment can be expressed as follows:  www.nature.com/scientificreports/ where p i is known, and σ R is the elastoplastic critical radial stress 21 . Thus, the radial stress components for the ith radius may be written as 21 According to the H-B yield criterion, the corresponding hoop stress can be expressed as: where H is defined in Eq. (4).
For the rock mass in the plastic zone, the elastic strain in terms of the stress considering the initial hydrostatic stress σ 0 , can be expressed as: Based on the plane strain condition, the elastic strain increment is calculated according to Hooke's law using the radial and tangential stress increments as follows 32 : Note that E (i) varies with the GSI and damage factor D and can be obtained using Eq. (12). where σ r = (σ (i) + σ (i−1) ) 2 . Therefore, the normalized radius ρ (i) can be written as: For the undamaged zone, the geometric equations in a polar coordinate system and the total strains of each annulus can be decomposed into elastic and plastic parts as follows 38 : Therefore, the compatibility Eq. (27) can be expressed as or (20)  where In the above equation, the dilation angle ϕ(η p ) varies with the increase in the softening parameter η p , and k (i−1) varies with the dilation angle. The deviatoric plastic shear strain η p (i) at the ith annulus can be updated as follows: Finally, the total radial and hoop strain at the ith annulus can be obtained as follows: Calculation of the stress-strain relationship in the damaged zone. In this study, we consider that the plastic region is divided into undamaged and damaged zones, as shown in Fig. 2, and assume that the stress and strain behaviors of the rock mass are continuous in the plastic zone. The solid red and blue lines represent the rock mass evolution functions in the undamaged and damaged zones, respectively. When the abovedescribed process is repeated m times, i.e., ρ (i=m) = ρ d . Subsequently, the surrounding rock is transformed from an undamaged zone to a damaged zone.
For the damage zone, the yield criterion can be expressed as follows 30 : As introduced above, because of the excavation-induced effects, the material parameters vary with the softening parameter and damage factor D in the damaged zone. For the H-B yield criterion, H d (σ r , η) can be expressed as: where m d ,s d , and α d are the strength parameters of the H-B criterion in the damaged zone, which can be calculated using GSI through Eqs. (42), (43), and (44).
The residual values of m, s, and α in the damaged zone can be rewritten as: www.nature.com/scientificreports/ where the subscripts "p" and "r" denote the peak and residual values of the material parameters in the damaged zone. When ρ (i=m) = ρ d , GSI d p = GSI (i=m) can be obtained, and GSI d r can be calculated using Eq. (11). As mentioned above, the stress and strain are continuous at the critical point in the damaged and undamaged zones. Under this condition, i = m represents the last annulus of the undamaged zone, and j(m ≤ j ≤ n) represents the number of annuli in the excavation damaged zone 25 . Thus, the radial stress for the jth annulus can be written as: and the hoop stress can be expressed as: The elastic strain increments in the radial and tangential directions in the damaged zone can be obtained using Eq. (24), and the Young's modulus is replaced with a new Young's modulus E d (j) . The governing equilibrium equations for the damaged zone are different from those for the undamaged zone given that the damage of the surrounding rock is considered. Therefore, the equilibrium equation and coordination equations under polar coordinates can be expressed as Hence, Eq. (50) can be approximated for the jth annulus as follows: where σ r(j) = (σ r(j) + σ r(j−1) ) 2 , and the normalized radius ρ d (j) can be expressed as follows: Because the strain is composed of elastic and plastic parts, Eq. (51) can be reformulated as: where ).
Similar to Eqs. (36) and (37), the radial plastic strain and the circumferential plastic strain in the damage zone can be obtained: . Subsequently, the plastic shear strain η can be expressed as: Hence, similar to Eq. (39), the total radial and tangential strains can be obtained: The normalized displacement can be determined using Eq. (56).
When the above calculation is iterated n times, the radial displacement of the tunnel surrounding rock is calculated as follows:

Verification and comparison
To prove the accuracy of the proposed method, which was implemented in MATLAB, the results obtained using the proposed solution are compared with those obtained by Mohammad Reza Zareifard 3 and Lee and Pietruszczak 21 .
The rock strength parameters corresponding to Zareifard's 3 study were taken as input data: b = 7 m, σ 0 = 27 MPa , p i = 5.14 MPa , σ c = 90 MPa , m i = 10,E 0 = 50 MPa , v = 0.25 (where E 0 is the Young's modulus of the intact rock). These data were used to derive the strength parameters of the rock masses in the different zones.

Model verification.
The accuracy of the algorithm can be verified by comparing its results with the existing results in literature without considering the disturbance damage of the surrounding rock.
The formation of the damaged zone around the tunnel is not considered, and hence, no excavation damagerelated formula was used for comparison. Figure 4a,b present the results of the two compared methods, which yield the same results as those obtained using the proposed algorithm, thus confirming the accuracy of the new user-coded algorithm.

Comparison with Lee and Pietruszczak's solutions.
To show the effectiveness of the proposed algorithm, the problem originally solved by Lee and Pietruszczak 21 was solved using the proposed method and then compared. The values of the damage factor D r were set to 0.2 and 0.4, and the results were compared with those obtained by Lee and Pietruszczak, who did not consider the damaged zone formation. Figure 5 shows the variation curves of the stress and radial displacement versus r/b after tunnel excavation. As shown, for D r ranging from 0 to 0.4, the radial displacement increases whereas the tangential stress decreases at the tunnel excavation surface.

Comparison with the results obtained by Zareifard and Fahimifar. A semi-numerical solution for
an axisymmetric circular tunnel excavated in a damaged rock mass is obtained using the algorithm proposed by Zareifard and Fahimifar 3 . The solution was proposed for solving circular tunnel excavation problems with a nonlinear H-B failure criterion. Figure 6 shows the distributions of the stresses and radial displacement calculated using both the approaches. Figure 6a shows that the distributions of σ θ and σ r are largely the same in the plastic residual zone, because the damage factor D r of the algorithm is equal to the damaged factor D of Zareifard and Fahimifar's algorithm for the residual zone. In Zareifard and Fahimifar's algorithm, the damage factor D is a constant value for the damaged zone, which can lead to sudden changes in the circumferential stress at the critical points in the elastoplastic zone. However, the authors proposed that the damaged factor D changes gradually from the elastoplastic critical point D p = 0 to the residual zone D r ; therefore, the σ θ distribution shows evident differences in the plastic zone. Meanwhile, the displacement calculated by Zareifard and Fahimifar is higher that calculated in this study.

Influence of new variables on strain-softening behavior
To illustrate the influence of the new variables on the results of the algorithm, the influence of each variable in the algorithm is studied in separate sections. In this case, the damage factor D and dimensionless damage radius ρ d are the main parameters characterizing the degree and extent of damage.
Effect of gradual variation damaged factor D. The parameter D r represents the maximum disturbance degree of the surrounding rock. Four common D r values, 0, 0.2, 0.4, and 0.6, were used to investigate the role of the damage factor D in the stress distribution and ground response curve (GRC). www.nature.com/scientificreports/ Figure 7a,b show the variation curves of the stress versus r/b in the elastoplastic zone for ρ d = 1 and ρ d = 0.9 . It is apparent that the circumferential stress increases gradually with r/b and reaches maximum at the elastoplastic boundary. As shown in Fig. 7, with the increase in the damaged factor D r , the circumferential stresses decrease on the surface of the tunnel opening, the inflection point between the plastic softening zone and the plastic residual zone is more evident, and the range of the plastic residual region increases significantly. As the rock mass damage occurs in the plastic zone (Fig. 7b), the circumferential stress shifts at the critical point between the damaged and undamaged zones. In addition, there is a remarkable drop in the distribution of σ θ in the damaged zone at higher damage factor values, which means that the increasing rate of the circumferential stress increases with the increase in the damage factor D r . Meanwhile, with the increase in D r , the slope of the radial stress curves decreases in the residual zone but increases in the softening area.
The evolution of the plastic radius R p versus the support pressure p i is plotted in Fig. 8, where a distinct difference in the evolution of the plastic radius can be observed for different damage factors D r . For different damage factors D r , the plastic radius R p calculated under the same ρ d condition increases with the decrease in the support pressure. Additionally, the plastic radius R p increased with the increase of damaged factor D r from 0 to 0.6. The GRC more intuitively presents the inverse relationship between the support pressure p i and the radial displacement at the tunnel excavation surface. Figure 9 depicts the GRC for different damaged factors D r . With the increase in the damage factor D r or the decrease in the support pressure, the radial displacement at the excavation boundary gradually increases. Furthermore, the effects of the damaged factor D r on the radial displacement are evident, particularly when the damage factor is D r = 0.6. Influence of damaged radius. As previously mentioned, the normalized damage radius ρ d represents the range of the surrounding rock disturbed by blasting or excavation. To investigate the influences of ρ d on the stress and GRC, we selected three different radii ρ d = 0.95, 0.9, and 0.85, and the damaged degree of the rock mass was set to 0.4. www.nature.com/scientificreports/ Figure 10 shows the distributions of the hoop stress σ θ and radial stress σ r with different normalized damage radii ρ d in the plastic region. As shown, the distribution of σ θ is largely the same in the residual zone, whereas it shows evident difference under different normalized damage radii ρ d in the softening zone. However, with the change in the damaged radius, the distribution of σ r shows no significant difference, and the tangential stress can be divided into two distinct stages: an undamaged zone stage and a damaged zone stage in the plastic-softening process. Thus, for ρ d ranging from 0.95 to 0.85, the transition points of the distribution of σ θ from the plasticsoftening zone to the residual zone and the undamaged zone to the damaged zone are evidently shifted to the left. However, the drop modulus of the circumferential stress increases with the decrease in the normalized damage radius in the damage zone.
The evolutions of the plastic radius R p and the radial displacement of the tunnel surface versus the support pressure are plotted in Figs. 11 and 12. Both these figures demonstrate a distinct difference in the evolution of the plastic radius and the GRC for different normalized damage radius ρ d . They reveal that the evolutions of the normalized damage radius have some influence on the plastic radius curve and GRC. With the decrease in ρ d , the plastic radius curve and GRC shift to the left, i.e., the radial displacements and plastic radius decrease under all internal support pressures.

Conclusions
In this study, the plastic zone comprises undamaged and damaged zones, where the damaged zone represents a gradual transition from a softening zone to a plastic residual zone. Based on the theory of progressive failure, a numerical solution for a circular tunnel excavated in a strain-softening rock mass was developed. In the finite extent damage zone, the reduced strength and stiffness of the surrounding rock with the gradual change in the damage factor D is considered for more accurate numerical results. The accuracy of the proposed algorithm, which takes advantage of two solution schemes (the iterative finite difference solution proposed by Lee and Pietruszczak and consideration of the damage zone as done by Zareifard 2020), was verified through some examples. The following conclusions can be drawn from the study: www.nature.com/scientificreports/ (1) By simplifying the algorithm as a special case, the results were found to be in good agreement with those of other methods. The algorithm could accurately reflect the strain-softening problem of a circular tunnel excavated in a rock mass. (2) While solving for the stress-strain states in the plastic zone, the variation in the Young's modulus proposed by Hoek, gradual variation in the damage factor D and normalized damage radius ρ d were considered, and the obtained results were compared with those of previous algorithms. Without considering the damaged factor during surrounding rock excavation, Lee's result found to be relatively conservative, and the damaged factor D in the damaged zone was constant, leading to an overestimation of the calculation result reported by Zareifard. (3) The comparative results showed that the variation in the damaged factor has distinct effects on the evolutions of σ θ , σ r , R p , and GRC in the plastic zone. The normalized damaged radius ρ d slightly affected the evolutions of σ θ , σ r , R p , and GRC.
In a nutshell: selecting an appropriate model that considers the damage radius and the variation in the damaged factor D is of great significance for analyzing the GRC of circular tunnels in rock masses exhibiting strain-softening behavior.   www.nature.com/scientificreports/